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ABSTRACT 


Combined  free  and  forced  convection  for  steady  fully  developed 
laminar  flow  in  long  horizontal  rectangular  channels  under  the  thermal 
boundary  conditions  of  constant  axial  temperature  gradient  and  peripherally 
uniform  wall  temperature  at  any  axial  position  is  approached  by  the  method 
of  successive-overrelaxation.  Two  vortices  (secondary  flow)  with  opposite 
sense  are  observed  in  a  cross  section  normal  to  the  main  flow  as  a  result 
of  the  buoyancy  effect.  The  convergence  of  the  numerical  solution  is 
ascertained.  Graphical  results  are  presented  for  velocity  and  temperature 
distributions,  streamlines,  isotherms,  w/Wq  vs.  ReRa,  fRe/(fRe)Q  vs.  ReRa 
and  Nu/NUq  vs.  ReRa  for  the  aspect  ratios  y  =  0.2,  0.5,  1,  2  and  5  and 
Pr  =  0.73.  For  a  square  channel,  the  heat  transfer  result  for  Pr  =  7.2 
is  also  presented. 
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NOMENCLATURE 


cross-sectional  area 
constant 

width  of  a  rectangular  channel 

constant,  2(1 /h^  +  1/J^) 

2 

constant,  1/n 

2 

constant,  l/£ 

height  of  a  rectangular  channel 

3  _ 

constant,  -  C-j  Dg  /4vy  =  Re/w 

axial  pressure  gradient,  9Pq/3Z 

axial  temperature  gradient,  9T/9Z 

equivalent  hydraulic  diameter,  4A/S 

Eckert  number  W^/C  C0D 

p  i  e 

friction  factor,  2Tw/(pW^),  or  a  function  of  X  and  Y 
functions  of  X  and  Y 

4  2 

Grashof  number,  £g  C2Dg  /v 
gravitational  acceleration 

dimensionless  grid  spacing  in  the  X  direction,  (r+l)/4M 

average  heat  transfer  coefficient,  -  q/T^ 

thermal  conductivity 

axial  direction  characteri stic  length 

dimensionless  grid  spacing  in  the  Y  direction,  (r+l)/(2yN) 


xi 


M 

N 

Nu 

n 


n 1 
Os 
P 

P1 


Pr 

Ra 

Re 

S 

T 


U,V,W 


u  ,v  ,w 


X  ,Y  ,Z 


x,y  ,z 


=  number  of  divisions  in  the  X  direction 

=  number  of  divisions  in  the  Y  direction 

=  Nusselt  number,  hD  /K 

e 

=  dimensionless  inward-drawn  normal  or  nth  iteration  for 
superscri pt 

=  dimensional  inward-drawn  normal 

=  Ostrach  number,  3gD  /C 

e  p 

=  pressure 

=  pressure  deviation  which  is  a  function  of  X  and  Y 

=  axial  pressure  distribution  measured  on  the  bottom  plate 

Y  =  0,  0  <_  x  <_  a  and  a  function  of  z  only 

=  pressure  at  wall,  Pn(z)  -  P  gY 

u  w 

=  Prandtl  number,  v/a 

4 

=  Rayleigh  number,  3gC 2D0  /va 
=  Reynolds  number,  Dg  W/v 
=  circumference  of  cross-section 
=  local  temperature 
=  mixed  mean  temperature  difference 
=  wall  temperature 

=  velocity  components  in  X,  Y  and  Z  directions 
=  dimensionless  velocity  components  in  X,  Y  and  Z  directions 
or  normalized  velocity  components  defined  in  (A.l) 

=  rectangular  coordinates 
=  dimensionless  rectangular  coordinates 


GREEK  LETTERS 


xn 


a 

e 

Y 

K 

e 

e 

0 


y 

y 

v 

P 


=  thermal  diffusivity 
=  coefficient  of  thermal  expansion 
=  aspect  ratio  of  a  rectangular  channel,  a/b 
=  vorticity  defined  in  (A. 3) 

=  a  prescribed  small  number 

=  small  quantities,  i  =  1,  2,  ...5  Defined  in  (A. 4) 

=  dimensionless  temperature  difference,  (T  -  T^J/fCgD  Pr  C) 
or  normalized  temperature  difference,  (T  -  ) /C2De 

=  viscosity 

=  a  value  defined  in  equation  (A. 8) 

=  kinematic  viscosity 
=  density 

=  mean  shearing  stress  at  wall 
=  dimensionless  stream  function 
=  relaxation  factor 
=  Laplacian  operator 
=  dimensionless  Laplacian  operator 


SUBSCRIPTS 


c 

i  J 
0 
W 


=  characteristic  quantity 

=  space  subscripts  of  grid  point  in  X  and  Y  directions 
=  condition  for  pure  forced  convection 
=  the  value  at  wal 1 


CHAPTER  I 


INTRODUCTION 


The  significant  effects  of  body  forces  on  forced  convective 
heat  transfer  in  internal  flows  have  been  studied  by  many  investigators 
in  recent  years.  Morton  [1]  analysed  the  buoyancy  effects  for  fully  de¬ 
veloped  laminar  forced  convection  in  uniformly  heated  horizontal  tubes 
at  low  Rayleigh  numbers  by  a  perturbation  method.  His  analysis  cannot 
be  readily  extended  to  the  cases  with  large  Rayleigh  numbers.  The  same 
problem  without  the  assumption  of  a  constant  axial  pressure  gradient  was 
discussed  by  del  Casal  and  Gill  [2]  using  a  perturbation  analysis.  A 
theoretical  investigation  of  natural  convection  effects  in  fully  de¬ 
veloped  horizontal  flow  between  infinite  parallel  plates  with  a  linear 
axial  temperature  distribution  was  presented  also  by  Gill  and  de  Casal  [3]. 
Mori  and  Futagami  [4]  reported  an  analytical  study  for  fully  developed 
forced  convective  heat  transfer  in  uniformly  heated  horizontal  tubes  with 
strong  secondary  flow  caused  by  buoyancy  forces  using  boundary  layer  approxi¬ 
mation.  Mori  and  his  co-workers  [5]  also  carried  out  an  experimental  study 
and  verified  the  result  of  theoretical  analysis.  The  influence  of  tube 
orientation  on  combined  free  and  forced  fully  developed  laminar  convec¬ 
tion  heat  transfer  was  discussed  by  Iqbal  and  Stachiewicz  [6]  using  a 
perturbation  approach  similar  to  Morton's  [1],  Recently,  an  experimental 


1 


2 


investigation  on  combined  free  and  forced  convection  in  a  horizontal  circu¬ 
lar  tube  was  reported  by  McComas  and  Eckert  [7]  with  several  earlier  experi¬ 
mental  investigations  quoted  in  the  references. 

It  is  noted  that  previous  theoretical  studies  on  combined  free 
and  forced  convection  in  horizontal  channel  flow  are  largely  confined  to 
circular  pipes.  In  contrast,  several  analytical  methods  are  available 
in  the  literature  for  combined  forced  and  free  convection  in  vertical 
channels  of  various  geometrical  shapes.  For  example,  combined  forced 
and  free  convection  for  fully  developed  laminar  flow  in  vertical  rec¬ 
tangular  channels  of  constant  axial  wall  temperature  gradient  with  uni¬ 
form  internal  heat  generation  was  treated  by  Han  [8]  using  double  Fourier 
series.  Tao  [9]  approached  the  same  problem  by  introducing  a  complex 
function  relating  to  the  velocity  and  temperature  fields. 

The  purpose  of  this  thesis  is  to  present  flow  and  heat  transfer 
results  for  fully  developed  laminar  forced  convection  with  helical  paths 
of  fluid  particles  due  to  the  buoyancy  effects  in  long  horizontal  rec¬ 
tangular  channels  having  aspect  ratios  0.2,  0.5,  1,  2  and  5  by  successive- 
overrelaxation  method  [10].  In  this  study,  the  wall  temperature  is  con¬ 
sidered  to  be  peripherally  invariant  and  a  linear  function  of  axial  position. 

Secondary  flow  for  convective  heat  transter  in  various  ducts 
is  also  known  to  occur  through  such  body  forces  as  centrifugal  and  coriolis 
forces.  Many  such  investigations  are  available  in  the  literature  [11-15]. 


CHAPTER  II 


THEORETICAL  ANALYSIS 


2 . 1  Governing  Equations  and  Boundary  Conditions 

Consider  the  steady  incompressible  fully  developed  laminar  flow 
in  a  uniformly  heated  horizontal  rectangular  channel  with  peripherally 
uniform  wall  temperature  at  any  axial  position.  The  motion  will  be  re¬ 
ferred  to  Cartesian  coordinate  (X,Y,Z)  as  shown  in  Fig.  1.  It  will  be 
assumed  that  the  effects  of  viscous  dissipation  (Os«l  and  Ec«l)  and 
of  compression  work  in  the  energy  equation  can  be  neglected,  that  vari¬ 
ations  in  density  due  to  temperature  changes  are  so  small  that  they  need 
to  be  taken  into  account  only  in  the  buoyancy  term,  and  that  the  kinematic 
viscosity  v  ,  and  the  thermal  conductivity,  K  ,  can  be  taken  as  constants. 
With  the  above  assumptions,  the  governing  equations  are  (see  APPENDIX  A): 
Continuity  Equation 


9U 

3X 


+ 


0) 


Momentum  Equations 


U  9U  + 
u  3X 


V 


au 

3Y 


1 

P 


9P 

3X 


+  v  V1  U 


(2) 


3 


4 


(4) 


It  is  noted  here  that  3W/3Z  =  0,  3U/3Z  =  0  and  3V/3Z  =  0  follow  from  the 
assumption  of  steady  fully  developed  flow. 

Energy  equation 


(5) 


where 


2  2 
3X^  3Y^ 


P  =  P 1 (X  ,Y)  +  Pw  ( Y ,Z)  =  P'  (X,Y)  +  PQ  (Z)  -  pw  g  Y 


The  buoyancy  force  in  equation  (3)  has  been  calculated  relative 


to  fluid  at  the  same  level  abjacent  to  the  channel  wall  and  the  remaining 
distribution  of  force  has  been  adsorbed  into  the  pressure,  P1.  For  steady 
convection  sufficiently  far  away  from  the  channel  opening  to  avoid  the 
entrance  length  effects,  the  temperature  throughout  the  fluid  increases 
linearly  along  the  axial  Z  axis.  Hence  the  secondary  flow  motion  driven 
by  the  buoyancy  forces  in  cross-sections  normal  to  the  channel  axis  is 
independent  of  axial  position. 


<f 
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The  boundary  conditions  are 

U=V=W=T  -  Ty  =  P '  =  0  at  channel  wall. 

2 . 2  Non-Dimensional  Transformations 

For  the  purposes  of  simplification  and  analysis,  the  above 
equations  (1  -  5)  may  be  reduced  to  non-dimensional  form  by  the  following 
transformations . 


X  =  De  x  ,  Y  =  De  y 
U  =  (v/De)u,  V  =  (v/De)v,  W  =  (vC/De)w, 

3PQ/3Z  =  C]  ,  3T/9Z  =  C2 
C  =  -  C^3/^,  T  -  Tw  =  (C2  De  Pr  C)0 

Re  =  WD  /v  =  C  w  (6) 

and  a  dimensionless  stream  function  \p  , 

u  =  8ip/3y,  v  =  -  dty/dx  (7) 


The  momentum  and  energy  equations  can  be  written  in  dimension- 


less  form  as  follows  after  eliminating  pressure  terms  between  equations 
(2)  and  (3),  and  applying  the  above  transformations. 
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3y  3x  v  ^  9x  9y 


(8) 


(9) 


(10) 


Although  the  partial  differential  equations  become  fourth  order,  the  5 
unknowns  U,  V,  W,  P1,  T  in  equation  (1-5)  are  now  reduced  to  3  unknowns 
tJj,  w  ,0  .  And  because  of  symmetry,  it  is  only  required  to  consider  half 
of  the  rectangular  region.  Consequently,  the  boundary  conditions  can  be 
rewritten  as: 


l^=-^7=w  =  0  =  Oat  channel  wall 


$  ■  § -  li  ■  i  ■ 0  at  x  =  a/2> 0  < y  <= b 


(id 


Solution  of  equations  (8  -  10)  satisfying  the  boundary  conditions 


(11)  is  a  matter  of  considerable  mathematical  difficulty.  Morton  [1] 
employed  a  perturbation  method  for  the  solution  of  the  corresponding  three 


7 


equations  and  the  associated  boundary  conditions  for  circular  pipes.  In 
this  study  for  rectangular  channels,  a  point  iterative  method  [10]  is 
used. 


CHAPTER  III 


FINITE-DIFFERENCE  APPROXIMATION 


3 . 1  Finite-Difference  Equations 

The  above  equations  (9)  and  (10)  can  be  regarded  as  having  the 

form 


^|-+^=F{x,y)  02) 

9x  9y 


which  is  called  Poisson's  (or  inhomogeneous  harmonic)  equation  and  its 
finite-difference  approximation  is  [10] 


02  .2 

f-’  -•  ~  p  2~  ^i  +  1  i  +  ^i  1  i^  +  2  2~  ^i  i+1  +  ^i  i-1^ 

i+‘  >J  i  i  jJ  2(h^+£^)  '  i  ' 


1  ,J  2(h2+£2) 


h202 

--/"T"  F(x,y) 


2(hc +£c) 


(13) 


The  inhomogeneous  term  F(x,y)  in  equation  (9)  can  be  approximated  as: 


F(x,y) 


dip  9w  9ip  9W  n 

9y  9x  9x  9y 


dip 


hi 


9y 


w.  ,1  •  —  W  •  -i  •  dip  •  • 
1+1  1-1  .  yi.J 


2h 


9x 


Wi  » j+1  ~  Wi  J-l 
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or 


F(x,y)  :  u.  . 

1  3  J 


wi>l,j  ~  wi-1,j  ,  v  Wi,itl'wi,j-1  1  (11) 

2h  vi  ,j  U  * 


Using  the  above  approximations,  the  finite-difference  equation  for  equation 
(9)  becomes. 


2  2 

wi,j :  [r  (wi+ij  +  +  r  (wij+i  +  "u-i) 


T-  ui  ,j  <wi+l,j  ■  w1-l,j)  '  Vvi,j  (wi,j+l  '  wi,j-l) 


+  2h2il2]/(h 2  +  J2) 


(15) 


Similarity,  equation  (10)  becomes 


ei,j  :  [2  (0i+l,j  +  6i-l,j)  +  2  (6i,j+l  + 


hH 


h2t, 


r Pr  uij  'Vi.j  ■  Vij)  -TPr  vij  (6ij+i 


2  2 

h  £  -i  w.2  ,  02n 

wij  — ]/(h  +  4  > 


(16) 


Equation  (8)  can  be  considered  as  inhomogeneous  biharmonic  equation  which 
has  the  form 


10 


V2V2  f  =  F'(x,y) 


or 


^  G  =  F(x,y) 


V  f  =  G 


Its  finite-difference  approximation  can  be  written  down  by  double  applica¬ 
tion  of  approximation  (13).  Then  equation  (8)  can  be  approximated  by 

,j  ~  ^2B1B2  ^i+l.j  +  'h-l  +  2B1B3  ^i.j+l  +  ^1  ,j-l ^ 


‘  B22  ^1+2 ,j  +  1('i-2,j*  "  B32  ^i,j+2  +  *i,j-2^ 


-  2B2B3  (<Fi+1j+1  +  Vi.j+i  +  Vl.j-l  + 


+  “k1  {B1  "  Vl.j)  +  B2  (*i+2,j  '  ^i-2 ,j 5 


+  B3  (Vl,j+l  +  Vl.j-l  ‘  Vl,j+1  ' 


v.  . 


+  “k1  {B1  '  ’h  ,j+l  ^  +  B3  ^i.j+2  "  ’h  ,j-2^ 
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+  B2  ^i+l.j+1  +  Vl.j+l  ‘  Vl.j-l  ‘ 


9  0  •  •  p  p  p 

+  Ra  C  ]/(B^  +  2B^  +  2B3  ) 


(17) 


where 


B]  =  2(B2  +  B3),  B2  =  1/h2  ,  B3  =  1 /l2 


30i,j/3x  ~  (0i-2,j  ‘  8ei-l ,j  +  86i+l,j  ‘  0i+2,j)/12h 

One  notes  that  the  finite-difference  approximations  for  the 
secondary  flow  velocity  components  u  and  v  are  as  follows: 


j  =  2, 


ui,2  ~  ^i,5  "  6^i,4  +  18^i,3  "  10^i,2^12A  O8) 


j  =  3,4,  ...,  N  -  1, 


ui  ,j  “  ,j-2  '  8^i  ,j-l  +  ^i.j+l  '  *i,j+2)/12,!'  ^19) 


j=N  u,  , 


^i,N-3  "  6^i  ,N-2  +  18<|,i,N-l  '  101,i  >N)/(-12i)  (20) 


where 


1  =  2,3, 


,  M  +  1 


r 
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i  =  2,  v2J  ~  (ip5  j  -  6^4jj  +  18*3J  -  10^2>j)/(-12h)  (21) 

i  =3,4,  . . . ,  M  +  1 , 


vi,j  '  ^i+2 ,j  "  ^l+l.j  +  8Vl,j  ‘  *i-2,j)/12h 


(22) 


where 


j  =  2,3, . .  N 

In  computation,  i  =  2,3,...,  M  +  1,  j  =  2,3,...,  N,  in  equations 
(15)  and  (16),  and  i  =  3,4,...,  M  ,  j  =  3,4,...,  N  -  1  in  equation  (17). 

The  boundary  conditions  (11),  take  the  following  finite-difference 

form: 

(a)  At  the  lower  horizontal  surface  0  <  X  <  a/2,  Y  =  0 


w.  =  6.  n  =  \K  ^  =  0 

1,1  1,1  M  ,1 


^i  ,2  ~  (3/4)  ^is3  ‘  (V3)  +  (1/16)  ip.  ^ 


(23) 


where 


i  =  1,2,.. 


M  +  1 
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(b)  At  the  upper  horizontal  surface  0  <  X  <  a/2,  Y  =  b 


wi,N+l  "  0i,N+l  ,N+1  '  0 


^i,N  ~  ^3/4^  ,N-1  "  ^1/3^  ,N-2  +  ^1/16^  ,N-3  ^24^ 


where 


i  =  1,2,...,  M  +  1 


(c)  At  the  vertical  surface  X  =  0,  0<Y<b 


w.  .=0.  .  =  ^ .  .=0 


^2, j  :  (3/4)  ^3J  -  0/3)  i|>4J  +  (1/16)  ip5jj 


(25) 


where 


j 


1,2 


.  ,  N  +  1 


(d)  At  the  vertical  middle  surface  X  =  a/2,  0  <  Y  <  b 


w 


M+2  ,j  5 


8M+2,j’  Vl.j  "  0 
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*M,j  =  '  *M+2,j>  =  '  *14+3,3 


(26) 


where 


j  =  2,3,. ,  N 

The  finite-difference  approximations  for  all  the  first  order 
partial  derivatives  in  this  study  are  obtained  by  trying  two,  three,  four 
and  five  grid  points  Taylor  series  expansion,  and  finally  five  grid  point 
approximations  were  found  to  be  reliable  and  satisfactory.  The  applica¬ 
tion  of  finite-difference  method  to  free  convection  problems  can  be  found 
in  references  [16  -  19]. 

3.2  Method  of  Solution 

(a)  Iterative  Procedure 

In  this  study,  the  successive-overrelaxation  method  discussed 
in  reference  [10]  was  employed  to  solve  a  set  of  finite-difference  equations 
(15  -  22)  and  the  associated  boundary  conditions  (23  -  26).  The  procedure 
is  di cussed  as  follows. 

1.  The  values  of  w.  .  at  each  node  point  are  obtained  by  solving 
equation  (15)  satisfying  boundary  condition  and  assuming 

u.  .  =  v.  .  =  0  for  all  i's  and  j's. 

1  jJ  1  )J 

2.  The  calculated  values  of  w.  .  and  u.  .  =  v.  .  =  0  are  used  in 

i  ,J  i  ,J  i  »J 

equation  (16)  to  compute  0.  .  at  all  interior  points  and  satisfy 

the  corresponding  boundary  conditions. 


«V  l 
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It  is  noted  that  the  above  two  steps  correspond  to  the  solution 

of  uniformly  heated  (or  cooled)  laminar  forced  convection  neglecting 

buoyancy  effect  due  to  the  variation  of  temperature.  The  values  of  w.  . 

1  sj 

and  0.  .  are  now  used  as  the  first  approximations  to  the  present  secondary 

•  jJ 

flow  problem. 

3.  Substituting  the  calculated  0.  .  into  equation  (17)  and  keeping 

'  5j 

u.  .  =  v.  •  =  0,  the  stream  function  ip.  .  at  each  node  point 

»  9J  I  Id  1  9d 

can  be  found. 

4.  Using  equations  (18  -  22)  the  values  for  secondary  flow  u.  . 

"I  jJ 

and  v.  •  can  be  obtained  for  all  the  node  points. 

5.  Using  the  values  of  u.  .  and  v.  .  acquired  in  step  4,  equations 

I  9j  1 

(15),  (16)  and  (17)  can  be  solved  again  for  w.  0.  .  and 

1  "I  3  J 

ip.  .  respectively.  The  new  values  of  u.  .  and  v.  .  can  be 

1  jJ  T  »J  1  >J 

calculated  from  the  latest  ip.  .  . 

6.  The  above  procedure  is  repeated  until  the  maximum  error  is  less 
than  a  prescribed  small  number  e  . 

The  details  of  solving  Poisson's  equation  and  the  inhomogeneous 
biharmonic  equation  by  a  point  iterative  method  are  to  be  discussed  in 
APPENDIX  B. 

(b)  Errors  and  Determination  of  Mesh  Size 

There  are  two  important  things  in  numerical  calculation.  One  is 
the  accuracy  of  the  results  and  the  other  is  the  computing  time.  In  the 
computation,  errors  can  come  from  the  following  sources: 


1 .  Round  off  error. 


l« 
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2.  Finite-difference  approximations  for  equations  (15  -  26). 

3.  The  convergence  of  each  Poisson's  and  inhomogeneous  bi harmonic 
equation  with  the  non-linear  terms. 

For  error  source  1,  because  of  double  precision  (16  Figures)  and  a  point 
iterative  method  used  in  this  study,  the  round  off  errors  should  not  be 
significant.  For  error  source  2,  this  inherent  error  can  be  reduced  by 
increasing  the  numbers  of  divisions  M  and  N.  A  grid  size  test  of  the 
convergence  of  the  numerical  method  was  made.  For  each  aspect  ratio  with 
Ra  C  =  0  and  10  ,  four  or  five  sets  of  M  and  N  (square  network  and  even 
numbers;  see  APPENDIX  B  for  details)  were  chosen.  Some  of  the  results  are 
shown  in  Table  7  and  Fig.  2.  The  number  of  divisions  adopted  are  also 
listed  in  Table  7.  For  error  source  3,  the  convergence  of  each  Poisson's 
and  inhomogeneous  bi harmonic  equation  including  the  non-linear  terms  is 
assured  by  continuing  the  computation  until  the  following  equation  is  satis¬ 
fied  . 


(n+1)  (n)  (n+1) 

f.  .  -  f.  .1  /If.  .  I 

l  ,j  l  ,j  1  max'  1  i  ,j  1  max 


(27) 


where 


e  =  5  x  10  for  each  Poisson's  and  inhomogeneous  bi harmonic 
equati on 

-5 

e  =  1  x  10  for  secondary  flow  u.  .  and  v.  .  . 

1  jJ  '  5J 

The  convergence  rate  is  faster  for  small  numbers  of  divisions  M 


3  T  - 


( 
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and  N.  In  order  to  ascertain  the  convergence  of  finite-difference  calcu¬ 
lations,  it  is  necessary  to  increase  the  values  for  M  and  N.  But  the 
larger  the  values  of  M  and  N,  the  slower  the  convergence  rate,  and  the 
greater  the  computing  time.  For  example  in  order  to  obtain  a  convergent 
solution  satisfying  condition  (27)  for  Ra  C  =  0  in  a  square  channel 
(y  =  1 )  with  M  =  16  and  N  =  32,  it  is  required  to  have  73  sweeps  for  the 
momentum  equation  and  56  sweeps  for  the  energy  equation.  But  if  M  =  10 
and  N  =  20,  it  only  takes  46  and  37  sweeps.  The  convergence  of  the  solu¬ 
tion  for  the  inhomogeneous  biharmonic  equation  is  slow.  Typically,  it 
takes  3  to  5  minutes  computing  time  for  a  complete  solution  for  a  number 
of  values  of  parameter  Ra  C  for  each  aspect  ratio  by  IBM  360  at  the 
University  of  Alberta  Computing  Center. 

The  expressions  F(x,y)  and  F'(x,y)  in  this  chapter  should  be 

F(x,y ’fx’ly)  and  F^x,y,fx  y2f  5ly  y2f^  ’  resPective1y-  These  lead  t0 
non-symmetric  matrices  for  linear  algebraic  equations. 


v,  2no  c  gsifqxn 
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CHAPTER  IV 


RESULTS  AND  DISCUSSION 

4.1  Flow  and  Heat  Transfer  Characteristics 

For  the  problem  under  consideration,  velocity  and  temperature 
distributions,  stream  function,  resistance  coefficient  and  Nusselt  number 
are  usually  of  interest.  Following  the  usual  definitions  the  expressions 
for  the  product  of  friction  factor  and  Reynolds  number  f»Re  and  the  Nusselt 
number  Nu  can  be  written  as  follows: 


f-Re 


W 


IP  w2 
2  w 


W  Do  p  2  t.  i  Da 
e  _  w  e 

y  y  W 


Nu 


(28) 


(29) 


where  and  h  can  be  found  directly  by  calculating  the  gradient  of 
velocity  and  temperature  on  the  boundary  r  respectively. 


,3W  n 

Tw  =  y 


w 


(30) 


-  K  l—l 

3n  W  W(T-T  )dA 

h  =  ’  TM  =  /A  W  dA 
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Substituting  equations  (30)  and  (31)  into  (28)  and  (29),  respectively, 
the  following  dimensionless  forms  can  be  obtained. 


f-Re 


2  (— ) 

^  ^3n;W 


w 


Nu 


36 

an 


w 


w 


(32) 


(33) 


Also  and  h  can  be  obtained  by  considering  overall  force  and  energy 
balances,  respectively. 


dZ 

31 

3Z 


A  r  A 
•  S  =  C1  s  = 


c 


C2  w 


D. 


4  T 


M 


(34) 


(35) 


The  results  are 


f*Re  =  2/w 
Mu  =  w2/4|w9| 

Equations  (32),  (33),  (36),  and  (37)  imply  that 


(— )  =  1 

v3n'w 


and 


30i 
3n  '  W 


(36) 

(37) 


=  w/4 


(38) 
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Some  of  these  numerical  values  are  listed  in  Table  8.  It 
can  be  seen  that  the  difference  between  these  two  methods  is  around  0.3% 
for  all  the  cases  considered.  The  values  given  in  Tables  1  -  6  are  com¬ 


puted  by  using  equation  (32)  and  (33).  The  above  average  quantities  such 
3w" 

3n  5  "w  3n 


as  w  ,  ,  w 0  and  are  obtained  by  using  Simpson's  rule. 


4.2  ReRa  Number  Effects 

To  illustrate  the  effect  of  the  product  of  Reynolds  and  Rayleigh 
number  on  the  axial  velocity  and  temperature  distributions,  the  profiles 
along  the  vertical  axis  X  =  a/2  are  plotted  for  the  aspect  ratio  y  =  1, 

Pr  =  0.73  in  Figs.  3  and  4.  It  is  seen  clearly  that  the  effect  of  the 
secondary  flow  is  to  shift  the  location  of  the  maximum  value  downward  and 
decrease  the  magnitude  of  the  maximum  value  as  the  parameter  ReRa  increases. 
The  velocity  and  temperature  distributions  along  the  horizontal  axis 
Y  =  b/2  for  y  =  1,  Pr  =  0.73  with  ReRa  as  a  parameter  are  shown  in  Fig.  5. 

As  expected,  the  magnitude  of  these  profiles  decreases  as  ReRa  increases. 

In  Fig.  6,streaml i nes  and  isotherms  are  plotted  for  a  square  channel  y  =  1 
with  Pr  =  0.73  and  ReRa  =  1.030  x  10^.  It  is  interesting  to  observe  that, 
due  to  the  higher  density  at  the  center  of  vortex,  the  location  of  the 
maximum  value  of  stream  function  ip  or  the  center  of  circulation  is  lower 
than  the  horizontal  central  axis,  Y  =  b/2.  One  also  notes  that  the  value 
of  stream  function  at  a  given  point  is  proportional  to  the  flow  rate  crossing 
any  path  joining  the  point  and  the  wall.  Although  not  shown  here,  the  value 


21 


of  the  stream  function  or  the  intensity  of  secondary  flow  increases  as  the 
value  of  the  parameter  ReRa  increases.  The  effect  of  buoyancy  force  on 
temperature  distribution  is  also  clearly  seen  in  the  distributions  of 
the  isotherms  on  the  right  hand  side  of  this  figure. 

4.3  Aspect  Ratio  Effects 

In  order  to  see  the  effect  of  aspect  ratio  on  flow  and  heat 
transfer  characteristics,  the  rectangular  channels  with  long  side  hori¬ 
zontal  will  be  considered  first.  The  velocity  and  temperature  distri¬ 
butions  along  the  vertical  axis  X  =  a/2  are  plotted  for  the  aspect  ratio 
y  =  2  and  Pr  =  0.73  with  ReRa  as  a  parameter  in  Figs.  7  and  8.  The 
velocity  and  temperature  profiles  along  the  horizontal  axis  Y  =  b/2  for 
y  =  2  and  Pr  =  0.73  are  shown  in  Fig.  9.  Nearly  the  same  trend  exists 
for  y  =  2  as  for  y  =  1.  Streamlines  and  isotherms  for  y  =  2,  Pr  =  0.73 

5 

and  ReRa  =  1.408  x  10  are  presented  in  Fig.  10.  For  this  case,  the  center 
of  circulation  is  seen  to  be  located  almost  along  the  horizontal  central 
axis  Y  =  b/2.  However,  the  location  of  minimum  temperature  (for  heating 
or  positive  ReRa)  is  still  found  below  the  horizontal  axis  Y  =  b/2.  The 
pattern  of  the  streamline  shows  that  the  intensity  of  the  secondary  flow 
near  the  vertical  wall  is  stronger  than  the  central  core  region. 

The  effect  of  aspect  ratio  can  further  be  seen  for  y  =  5.  The 
velocity  and  temperature  profiles  along  the  vertical  central  axis  X  =  a/2 
are  plotted  in  Figs.  11  and  12,  respectively ,  with  Pr  =  0.73  and  ReRa  =  0 

5 

and  2.984  x  10  .  The  velocity  and  temperature  profiles  along  the  hori- 
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zontal  central  axis  Y  =  b/2  are  shown  in  Fig.  13.  These  figures  show 
that  the  secondary  flow  has  less  effect  on  the  velocity  and  temperature 
distributions  in  the  central  region  for  higher  aspect  ratio  y  .  This 
phenomenon  can  be  seen  also  from  Fig.  14  where  the  streamlines  and  iso¬ 
therms  are  shown  for  y  =  5.  It  is  seen  that  the  secondary  flow  is  weak 
in  the  central  region  and  rather  intense  near  the  short  vertical  wall. 

The  center  of  circulation  for  this  case  is  almost  along  the  central 
horizontal  axis  Y  =  b/2.  The  streamlines  for  the  rectangular  channels 
with  y  =  1,  2  and  5  show  clearly  that  as  the  aspect  ratio  y  increases, 
the  location  of  the  center  of  circulation  tends  to  more  toward  the  vertical 
short  wall.  As  a  result  of  this,  the  secondary  flow  is  rather  weak  in 
central  region  and  intense  near  the  vertical  wall.  It  is  noted  that  for 
y  =  2  and  5,  the  center  of  circulation  is  almost  along  the  horizontal 
central  axis  whereas  for  y  =  >  the  center  of  circulation  is  well  below 

5 

the  central  horizontal  axis  for  ReRa  =  1.030  x  10  . 

The  effect  of  aspect  ratio  when  the  long  side  is  vertical  will 
be  considered  next.  The  velocity  and  temperature  distributions  along 
central  axes  X  =  a/2  and  Y  =  b/2  are  shown  in  Figs.  15,  16  and  17  for  a 
rectangular  channel  with  aspect  ratio  y  =  0.5,  Pr  =  0.73  with  ReRa  as  a 
parameter.  Comparing  with  the  case  y  =  2,  one  notes  that  in  general  the 
location  of  the  maximum  value  for  velocity  and  temperature  profiles  moves 
closer  to  the  bottom  horizontal  wall.  The  streamline  and  isotherms  for 

5 

the  aspect  ratio  y  =  0.5,  Pr  =  0.73  and  ReRa  =  1.381  x  10  are  shown  in 
Fig.  18.  One  sees  that  the  center  of  circulation  is  well  below  the  hori- 
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zontal  central  axis  and  the  secondary  flow  motion  is  weak  in  the  upper 
region  and  rather  intense  near  the  bottom  horizontal  wall.  Because  of 
rather  intense  secondary  motion,  large  temperature  gradients  exist  near 
the  bottom  horizontal  wall. 

The  graphical  results  for  velocity  and  temperature  distributions 
for  y  =  0.2,  Pr  =  0.73,  ReRa  =  0,  2.535  x  10^  are  shown  in  Figs.  19,  20 
and  21.  For  this  case  the  location  of  the  maximum  magnitude  for  both 
velocity  and  temperature  is  seen  to  be  quite  close  to  the  bottom  hori¬ 
zontal  wall.  The  streamlines  and  isotherms  for  y  =  0.2,  Pr  =  0.73  and 

5 

ReRa  =  2.535  x  10  are  shown  in  Fig.  22.  Once  again  one  sees  that  stronger 
secondary  motion  exists  near  the  narrow  bottom  horizontal  wall.  Corres¬ 
pondingly,  larger  temperature  gradients  exist  there.  Comparing  with  the 
case  y  =  0.5,  one  notes  that  the  center  of  circulation  and  location  of 
the  minimum  magnitude  of  temperature  move  further  downward. 

4.4  Results  of  w/w~q,  fRe/lfRe)^,  Nu/Nuq  and  Pr  Number  Effect 

It  may  be  of  some  interest  to  compare  the  relations  between 
w/wQ  and  ReRa  for  various  aspect  ratios  considered.  The  results  are  shown 
in  Fig.  23  for  Pr  =  0.73.  The  effect  of  secondary  flow  on  the  average 
axial  velocity  is  greatest  for  y  =  1  and  decrease  as  y  ■>  0  or  »  .  A  plot 
of  the  ratio  of  the  product  of  friction  factor  and  Reynolds  number  with 
and  without  buoyancy  effect  fRe/(fRe)Q  versus  the  parameter  ReRa  is  pre¬ 
sented  in  Fig.  24  for  Pr  =  0.73  with  aspect  ratio  y  as  a  parameter.  A 
similar  plot  for  the  ratio  of  Nusselt  number  with  and  without  buoyancy 


24 


effect  Nu/Nuq  is  shown  in  Fig.  25.  In  order  to  see  the  effect  of  Prandtl 
number  on  heat  transfer,  the  result  for  y  =  1,  and  Pr  =  7.2  is  also  shown 
in  the  same  figure.  Since  the  effects  of  Prandtl  number  on  the  ratios 
w/wq  and  fRe/(fRe)g  are  not  very  significant,  the  results  for  the  case 
Pr  =  7.2  are  not  shown  in  Fig.  23  and  24  (see  Table  2  and  equations 
(8  -  10)  also) . 

4.5  Conclusions 

1.  The  method  of  successive-overrelaxation  can  be  applied  to 

combined  free  and  forced  laminar  convection  in  long  horizontal 

rectangular  channels  with  various  aspect  ratios.  Unfortunately , 

this  method  diverges  for  the  values  of  the  parameter  ReRa 

beyond  those  shown  in  Figs.  23  -  25  and  Tables  1  -  6.  The 

exact  reasons  for  the  divergence  of  the  solution  were  unknown 

to  the  author.  However,  it  was  found  that  the  numerical  values 
2  2  2  2 

of  the  terms  V  V  ip,  V  w  and  V  0  approach  zero  as  the  values  of 

ReRa  becomes  larger  than  indicated.  This  effect  can  also  be 

seen  in  normalization  analysis  (see  APPENDIX  A).  These  terms 

1  /2 

are  all  of  the  order  of  1/Gr  .  Despite  this  difficulty,  the 
method  yields  a  satisfactory  solution  for  a  range  of  ReRa  far 
exceeding  that  of  Morton's  perturbation  method  [1]  for  laminar 
convection  in  uniformly  heated  horizontal  pipes  in  error  by 
about  10%  for  RaC  =  3.84  x  10^. 
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2.  As  can  be  expected,  the  resistance  coefficient  and  the  Nusselt 
number  increase  as  the  value  of  ReRa,  representing  buoyancy 
effect,  increases.  The  effect  of  aspect  ratio  on  flow  and  heat 
transfer  characteristics  is  clearly  seen  in  Figs.  24  and  25. 

For  the  aspect  ratio  y  =  1,  the  effect  of  Prandtl  number  is 
considerable.  As  similar  effect  can  be  expected  for  other  aspect 
ratios . 

3.  The  details  of  secondary  motion  caused  by  buoyancy  forces  are 
clearly  seen  from  the  pattern  of  streamlines  from  which  it  is 
evident  that  there  are  two  vortices  with  opposite  sense  formed. 

It  is  also  noted  that  for  a  square  channel  the  center  of  circu¬ 
lation  is  well  below  the  central  horizontal  axis.  One  can 
visualize  the  secondary  flow  motion  by  noting  that  the  flow 
rate  through  any  line  connecting  the  center  of  circulation  and 
a  point  on  the  wall  is  constant  for  a  given  value  of  the  para¬ 
meter  ReRa.  For  a  rectangular  channel,  with  long  side  horizontal, 
as  the  aspect  ratio  y  increases,  the  intensity  of  secondary 
motion  near  the  short  vertical  wall  increases.  For  a  rectangular 
channel  with  long  side  vertical,  as  the  aspect  ratio  decreases, 
the  center  of  circulation  moves  toward  the  bottom  short  hori¬ 
zontal  wall.  Consequently,  intense  secondary  motion  exists  near 
the  bottom  wal 1 . 

4.  The  method  used  here  has  also  been  applied  to  a  circular  pipe 
as  Morton  did  in  [1].  But  due  to  the  slow  convergence  of  the 
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inhomogeneous  bi harmonic  equation  in  polar  coordinates  and  the 
difficulty  with  boundary  conditions  for  the  vorticity  function, 
i t  did  not  work  wel 1 . 

5.  The  perturbation  method  used  by  Morton  may  be  applied  here  if 
numerical  integration  is  used  instead  of  an  analytical  method 
of  solution.  This  possibility  is  discussed  in  APPENDIX  C. 
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APPENDIX  A 


NORMALIZATION  OF  EQUATIONS 

To  normalize  the  governing  equations,  first  introduce  the  fol¬ 
lowing  dimensionless  quantities, 

X  =  D0x  ,  Y  =  Dey  ,  Z  =  Lz,  U  =  Vcu,  V  =  Vcv 
W  =  Ww  ,  PQ  =  PcpQ  ,  T  -  Tw  =  C2De0  (A. 1 ) 

The  pressure  terms  can  be  eliminated  between  the  momentum 
equations  in  X  and  Y  directions.  Noting  that  the  driving  forces  (axial 
pressure  gradient  and  buoyancy  forces)  and  the  inertia  terms  can  be  con¬ 
sidered  as  the  same  order  of  magnitude  and  defining  characteristic 
quanti ti es , 


Vc  =  (eg  C2De)1/2  =  Gr1/2  v/De  , 

P  =  v  p  W  Gr1/2L/De  =  p  W2  (Gr1/2/Re) (L/D  )  (A. 2) 

then  the  governing  equations  can  be  written  as  follows, 
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3u/3x  +  3v/3y  +  e-j  3w/3z  =  0 
E,  =  3u/3y  -  3v/3x 

?  i  /?  ?  ? 

u  3£/3x  +  v  3£/3y  +  e-j  w  3£/3z  =  V  ?/Gr  7  +  e2  3  £/3z  -  30/3x 

21/2  22 

u  3w/3x  +  v  3w/3y  +  e-|  w  3w/3z  =  -  3Pg/3z  +  V  w/Gr  7  +  e2  8  w/3z 

u  30/3x  +  v  30/3y  +  e]  w  30/3z  =  V20/(Pr  Gr1/2)  +  e2  320/3z2  +  $ 

0  =  e3  [2{ ( 3u/3x) 2  +  (3v/3y)2}  +  (3v/3x  +  3u/3y)2 

+  (De/L)2  { ( 3v/3z) 2  +  (3u/3z)2}] 

+  e4  [2(De/L ) 2  ( 3w/3z) 2  +  (3w/3y)2  +  (3w/3x)2] 

+  ec  [2(D  /L)2  { ( 3w/3y )  (3v/3z)  +  (3u/3z)  ( 3w/3x) }]  (A. 3) 

o  e 

where 

e,  =  (Re/Gr1/2)  (De/L)  ,  e2  =  { Dg/ L ) 2/Gr 1 /2 

e,  =  Os/Gr1/2  ,  e4  =  Ec/Gr1/2  ,  e5  =  Ec/Re  (A. 4) 


32 


Since  the  flow  is  fully  developed,  De/L  <  <  1  and  viscous  dissipations 
are  negligible  i.e.  Os  <  <  1  and  Ec  <  <  1.  Thus,  all  the  terms  with 
coefficients  e. (i  =  1,2,... 5)  can  be  neglected,  and  the  above  equations 
(A. 3)  can  be  written  as, 


3u/3x  +  3v/3y  =  0 


£  =  3u/3y  -  3v/3x 

u  3£/3x  +  v  3£/3y  =  V2£/Gr1/2  -  30/3x 
u  3w/3x  +  v  3w/3y  =  -  3Pq/3z  +  V2w/Gr^2 

u  36/3x  +  v  30/ 3y  =  V20/(Pr  Gr1/2)  (A. 5) 

It  is  noted  that  if  the  value  of  Gr  number  is  less  than  the 
order  of  102,  then  the  terms  V2£/Gr^2  ,  V2co/Gr^2  and  V20/(Pr  Gr^2^ 
can  not  be  neglected.  In  contrast,  Mori  et  al  [4]  neglected  these  terms 
in  their  analysis  for  higher  value  of  Gr  number. 


APPENDIX  B 


METHOD  OF  SOLVING  EQUATIONS  V2f  =  F  AND  V4f  =  F' 

It  is  desirable  to  discuss  a  procedure  for  solving  the  Poisson's 

2  4 

equation  V  f  =  F(x,y)  and  inhomogeneous  biharmonic  equation  V  f  =  F'(x,y) 

with  the  associated  boundary  conditions  by  a  point  iterative  method  which 

was  mentioned  in  3.2.  The  finite-difference  approximation  of  Poisson's 

equation  (13)  can  be  rewritten  as 


i 


Al(fi+l,j  +  +  A2(fi,j+1  +  T'J-l*  '  A3  Fi,j  (BJ) 


where 


Ai  ■ 


i 


2(h2+£2) 


A2  = 


2  ( h2+il2 ) 


A3 


2  2 

hV 


2(h2+£2) 


At  the  beginning,  prescribe  the  value  of  boundary  points  and  let  f.  .  =  0 
for  all  i's  and  j's  in  the  region.  The  calculation  is  executed  from  i  =  2, 
j  =  2  along  the  column  i  =  2.  For  example,  the  value  of  f2  2  1S  computed 

by 


f2 ,2  A1  +  fl  ,2^  +  A2  +  f2 ,1 ^  '  A3  F2 ,2 


(B .  2) 
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The  newly  obtained  value  2  1S  usec*  in  the  computation  of  f2  3  by 

f 2 ,3  =  Al^°  +  fl  ,3^  +  M°  +  f2,2^  "  A3F2 ,3  ^B*3^ 

This  process  is  continued  up  to  j  =  N  and  then  similar  process  is  re¬ 
ported  for  a  new  column  i  =  3  and  the  succeeding  columns  until  the  last 
column  i  =  M  +  1  is  reached.  After  completing  one  whole  sweep,  the  calcu¬ 
lation  is  back  to  i  =  2,  j  =2  again.  For  example, 


f2  2  =  Mf3  2  +  fl  2)  +  Mf2  3  +  f21^"A3F22  (B-4) 


where  the  values  of  f^  2  anc*  ^2  3  come  ^rom  t*ie  resu^s  of  the  first 

sweep.  The  above  process  is  repeated  until  equation  (27)  is  satisfied. 

Then  a  numerical  solution  f.  .  which  satisfy  the  Poisson's  equation  and 

the  associated  boundary  condition  is  obtained. 

4  /  v 

The  method  of  solution  for  V  f  =  F'(x,y)  is  similar  to  the  above 
except  that  the  computation  is  started  from  i  =  3,  j  =3  and  the  values 
of  boundary  and  adjacent  points  are  obtained  by  two  boundary  conditions. 

In  carrying  out  the  above  calculations,  the  following  formula 
for  successive-overrelaxation  is  used  [10]: 


f (n+1 ) 
i  ,j 


+ 


f (n+1 ) 
Ti-1  ,j 


)  +  A2(f 


(n) 
i ,  j+1 


+  f 


(n+1) 
i  »j-l 


1 


(B.5) 


The  superscript  (n)  denotes  that  the  values  obtained  by  the  previous 
sweep.  The  relaxation  factor  w  is  calculated  by  the  following  equation 
DO], 
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u>  =  1  +  [y/{l  +  (1  -  jO  '  )] 


A 1/2,12 


(B.6) 


where 


(B.  7) 


a/2  =  Mh  ,  b  =  Nh  and  h  =  mesh  size. 

It  is  noted  that  the  square  network  h  =  £  was  used  for  all  the  cases 
considered  in  this  study  and  even  numbers  of  M  and  N  were  taken  for  the 
purpose  of  applying  Simpson's  numerical  integration.  The  values  for 
the  relaxation  factor  oj  calculated  by  equation  (B.6)  are  optimum  only 
for  essentially  self-adjoint  equations. 


APPENDIX  C 


PERTURBATION  METHOD 

Solution  of  equation  (8  -  10)  and  boundary  condition  (11)  can 
also  be  obtained  by  expanding  ip  ,  w,  and  0  as  power  series  in  the  product 
of  Rayleigh  number  and  C,  provided  that  this  is  numerically  small. 

ip  =  (RaC)  i|j-|  +  (RaC)2  ip^  +  ... 

w  =  w  +  ( RaC  )w  1  +  ( RaC ) 2  w  2  +  ... 

9  =  eQ  +  ( RaC ) 0-j  +  (RaC)2  02  +  ...  (C.l) 

The  leading  term  of  ip  must  vanish  because  there  is  no  circulation 
when  buoyancy  effect  is  neglected  or  RaC  =  0.  After  substituting  the 
above  power  series  into  equation  (8  -  10),  ip.  ,  w^  and  e.  are  obtained 
by  equating  coefficients  of  the  same  powers  of  RaC.  The  coefficients  of 
(RaC  ft  are 


V2wQ  +4=0 

y20o  +  wQ  =  0  (C .  2) 
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The  coefficients  of  (RaC)  are 


=  0 


Sip-j  3Wq  3^i  3Wg  2 
37”  3x“  "  3x“  3y~  =  V  W1 


Pr  ( 


3^1  30q  dip 


1  360 


3y  3x  3x  3y 


)  =  V^Q-j  -  w-j 


(C .  3) 


and  so  on. 

It  is  quite  tedious  to  obtain  an  analytic  solution  of  the 
above  equations  in  a  rectangular  channel.  Thus  any  numerical  methods 
can  be  used  to  solve  them.  Due  to  the  existence  of  errors  in  numerical 
solutions  of  (C.2),  (C.3)  »and,  recalling  the  product  of  (RaC)n  and 
(or  wn»  6  )  in  Series  (C.l),  the  error  of  the  solution  ip  (or  w,  0) 
become  uncontrollable  when  the  value  of  RaC  becomes  larger.  The  rather 
limited  ranges  of  RaC  and  the  unknown  errors  are  the  reasons  that  this 
method  were  not  tried  here. 


APPFNDIX  D 


FORTRAN  PROGRAMMING 


DOURLF  PRECISION  PF  T  A  ,  FRR  » FM  »FN  *  FI »  L  9  P 1 » P2  *  R?  *  L?  >  HR  * 

1  W  (  9  -a  ,  1  0  H  )  »T  (-2^*10?)  .  F  (  7.  2  ,  ]  n  0  )  ,  U  (  p  p  *  1  02')  »V(  22 , 1  02  )  » 

?GAWMA*PR  ,RAr»PJ  *PFF  9  W  m  ,  T  M  A '  >  FR  F  •>  T  F  *NU*FR1  ?  V  T  E  S  T  * 

pwa  ,WN  ,WD  ?  WFR  ,  TA  *  TN  %  TO  *HA  *.TFR  *RCT  *  SA  »  $A  ,SFR  ♦  VN  •  UN  *UD  9  \/D  » 

4SH  (  PR  9 1  OP  )  ,  TG  (  PP  9 1  OP  )  ,SIJM  *  SUT  *SUM3  *  5UT1  *  SHM  ?  TGM  *  HP  *VFR  ,  S'  , 
PCI  > f  2  » C P  » C4  9 C 5  > C 6  *  4  »  X 1  >  I  ?  *XP  »X4»X5  >X6.*X7.»X8  *X9  *X1  -  1 

I N TEGFR  T IM1 *T  I MS  ,WK * TK * SK 

m  IF  IMF  NO.  OF  PIVIFIONF  IN  THE  X  DIRECTION  BETWEEN  <"»  ,A/2 

N  IF  THF  NO.  OF  D  T  V  T  f I o  N  S  IN  THF  Y  DIRECTION  R  F  T  W  F  E  N  0  »R 

PETA  TF  THF  RELAXATION  FACTOR 

ERR  IF  THE  PRESCRIBED  ALLOWABLE  ERROR 

TPM  AND  TI-VP  ARF  THF  PRFFCRIpED  NO.  OF  FWFFPp 

R F  A  n ( F , 4  p  )  M  >M jB^TA  >CRR  .TIM!  9 T I M ?  *  T I  V P 

Ml -M+ 1 

M?=M+? 

M  “2  =  M  +  P 
M  I  =  M  -  1 
Ml  =  N  +  1 

M  I  =  (\l  -  1 

N  J  -  N  —  ? 

NK =  N  —  3 

fm  =  drl.f  (  float  ( m ) ) 

FM  =  DRLE ( FLOAT ( N  )  ) 

U  v  w  ARE  THE  DIMFNFIONLFSS  VpLOCITY  COMP.  TN  X  Y  Z  DIRE. 

T  IF  THp  n IMENF tonlffs  tfmp. 

F  IF  THF  FTREAMFliNCTlON 

I  J  ARE  THE  SPACE  SUBSCRIPTS  OF  GRID  PT.  IN  X  Y  DIRE* 

DO  i  1  =  1  ,  M  P 
DO  1  J=1 *N1 
m (  I  ,  J ) =DDO 
t (  1  ,  J ) =nnn 
p (  I  •  J ) =  PDA 
U (  I  ,  J)=DDF 
I  V (  I  ,  J ) =^DF 
S 5 F  CONTINUE 

K  Ip  THF  NO.  OF  SWEEPS  OF  IJ  V 

GAMMA  IS  THF  asPFft  RATIO  OF  A  RFC T ANGULAR  CHANNEL 
PR  IF  THF  PRANDTL.  NO. 

RAT  IF  THE  PRODUCT  OF  RA  AND  C 

H  L  ARF  THE  01 MENSTONLFSS  GRID  SPACING  IN  THF  x  Y  DIRE. 

K  =  1 

RFAD(S  9  46  )  G  A  M  Vi  A  PR  ,  R  A  f 
H=  (  GAMMA+1  Da  )  /  (  4D0*F.M  ) 

|_  =  (  GAMMA  +  1  Dp  )  /  (  ?UP*GAMMA#FN  ^ 

L2=L*L 


H2=l  !#H 

r  '  =  ]  j  /■'? 

r'^  =  1  DL/L2 

n  i  =pnr*  { !;,7+r  7 ) 

r  4  =  o  i  *  7  +2  DO  *  p  ^  -* -ij  p  +  p  p  fl#  7  7  #  --  7 
Ci  =  7n-r-*L2/  (  H2  +  L.7  ) 

0  7  =  ^  D  - 1 4  H  ?  / ( H2  +  L2 ) 

r  ^  =  -  n  - 1  *  h  c  1 

C4  =  £'P-l*  L#C? 

CF^HP^Cl 
C6=H5#4D0 
X1=2D(.  *  B 1  #  HS  ?  /  R  4 
X2=7D0*R1#33/Q.4 
X^=°.?*^?/D4 
X4  =  n  3**2  /R4 
X  5  =  2.  D  0*  B  2"$  B  3 / B4 
X6  =  r>D- 1  *  1  /  (  H*jB4  ) 

X  7  = t;  D  —  1  »  n  7  /  ( i  :  !!‘n  4  ) 

X  0  =  S D- 1 "  '°/tH*P4) 

X9»FD*-1*B1 / ( L*B4 ) 

XI p  =  5D-l *P?/ ( L  * P  4 ) 

XI  i  =5D-1  *B3/35t^B4  ) 

WRITE  (6  $4° )  ' '  •  -  "  ’  *  '  *  » PR  * 

t  is  THF  NO ♦  :  SWEEPS  Of  FHf  "  :  Tt :  '  EON. 

1  T  1  I  l  1 1 ! 1  .  v  I  CITY  PROFILE 


'K  = 

1 

n  T  = 

ODD 

Rfr 

=  7  D  7 

DO 

4  1  = 

2  *  ‘ '  1 

11  = 

1  +  1 

I  J  = 

1-1 

DO 

4  J  = 

2  9  N 

Jl  = 

J  + 1 

JI  = 

J-l 

’  A  = 

^  1  *  ( 

W  (  1 1  » J  )  +  (  I  I  >  J  )  )  +  C2  #  <  W  (  T  *  J 1  )  +W  (  T  *  J  I 

\  (  v  (  n  ,  J  )  -!■;  (  I  I  ,  J  )  )  -C4*V  (  I  » J  )  *  (  W  (  I  ,  J  1  )  (  I  ,  J  I  )  )  +CS 

’,•71  =tDFT  A  *W  A  -  (  c  ~  j  a  -1  p  n  )  *  Y  (  I.»  J  ) 

I F ( I . LT . y .OR . I • GT . M )  CO  TO  5 
{  w  7  ,  j  )  =  wfv 

r  ’..I fp  =  p  A  07  ( -v'  (  I  *  J  )  ) 

PT  j  c,  7  per  rs/  A  x  •  D  I  F  F  r  R  E  f\|C  P  °  ~  jucrr  \|  JH^  VALUE  J I  1ST 
c  OPTATNFD  AND  PREVIOUS  VALUE 

t  =nv  a  xi  (pi,  •’P  ) 

C;  REF  I S  T 1  IF  VAX*  AQ  SOLU  TE  VALIJF 

7EF  =  DUAX1  (  °.pr'  ) 

4  u  ( i ♦  j  )  ='*'n: 

F  R  =  D  I  /  orF 

IF  (u'FR.LE.ERR.OR.  W  .Or.Tl  '1  )  00  TO  6 

v  kr  =  iv  v  + 1 
GO  TO  7 


6  C^MTIMUF 

ti<  T  S  THF  MO.  OF  S'a,F F  P  S  OF  THF  FNFRGY  FON* 
lie  TMG  THC  OR  TAIN  F  D  TO  F  I  NR  tFMP.  PROFILE 
tf  -  i 

11  ^  T  =OnO 

nFFrinn 
OO  P.  jM  i 

T  1  =  T  +  i 
T  T  =  T  -1 
00  R  j=? ,N 
Jl =  J  +  1. 

J I  =  J  - 1 

T  A  =  r  i  *  (  T  (  T  i  ,  J  )  +T  (  T  T  ,  J  )  ) +07* ( T (  I  9  J]  )  +  T  (  I  *  J  I  )  )  -PR*  (  oo*[i  (  T  ,  J  )  * 
1  ( T  (  t i  ,  J  )  -T ( t  T  ,  J  )  ) +  r 4 * V ( I > J ) * ( T ( I  ,  J  l  )-T (  I  »  J  I  )  )  ) -  V ( I  * J) *CS 
TISJ  =  nFTA*TA-(PFTA-1O'0)*T(  T  «  J  ) 

IF (  I • L T • M . OR • I •GT.M)  00  TO  9 
T  (  m  ,  J  )  =  T  N 

o  j n  =  n  a D  p ( TM-T (  I  » J  )  ) 

O  t  =om a xi  { D I  * TD ) 

HA=0APS(TN) 

RFF  =  DY.A  X  1  (  REF  ,  HA  ) 
p  t  (  T  .  J  )  =  T  N 
TF7  =  n i  /p^f 

IF ( TFR.LF.FRP.OP.TK .0F.TTM1 )  00  TO  10 

TF  =  tf  +  ] 

On  TO  1 1 
1  n  r  O  m  t  T  m  i 1 F 

PF  TP  THF  MO.  OF  PWFFPP  OF  THF  VO  R  T IT  T  T  Y ( F  T  R  F  A M  F !  J  N C  T I O  M ) F  QN 
iipjMO  T HO  ORTAIMFO  T  TO  OOMPUTF  STRFAMFlJMOT  TOM 
PK  =  1 

IF  PI  =  '••'00 
p  F  c  =  o  O  0 

DO  12  1=3?  M 
I  i  =  t  +1 
I  2  =  I  +  2 
IT  =  T-1 
I  J  =  I  —  7 

no  1  7  J  =  7  ,  m  j 
Jl =J+1 
J  2  =  J  +  2 
J  T  =  J  - 1 
J  J  =  J  —  2 

ROT  =  R  AC*  (  T  (  T  J  >  J  )  -PPP*T  (  I  T  » J  )  +RDO*T  (  T 1  » J. )  -T  (  T.?  *  J  )  )./  (  1  7D0*H  ) 

s A  =  xi * ( s (  II  ,  J ) +F (  IT  , J )  )  +  X  7  * ( S (  T  , Jl  ) +S {  T  9 Jl  >  )-X7* 

1  (  S  (  T  7  9 J ) +S (  T J , J  )  ) 

1  —  X  4*  {  S  (  I  9  J7  )  +s  (  I  9  J  J  )  )  -X  5*  (  S  (  II  ,  J  1  I  +S  (  I  f » J 1  )  +.s  (  I  I  9  J  I)  + 

1  f  (  T 1  ?  J 1  )  )  + 

7U(  I  9 J 1* ( X6* ( S ( I  I  * J )-$ ( 1 1  9 J )  ) +  X7* ( S ( 12  9 J )-S (  I  J  * J  >  ) 
p  +  X  R* ( S (  II  9  J 1  ) +  S (  I  1  9  J  T  ) - S ( I  I  . J 1  ) - S (  I  I  ,  J  I  )  )  ) 

4  +  V  (  T  ,  J  )  *  LXQ*  (  S  (  I  9  J  T  )  -S  IT  9  Jl  )  )  +X  1 1*  (  S  (  I  9  J'2  )  -S  (  I  »  JJ  }  ) 
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■■  i  ■■  1  '  '  ■;  :  " 1  . .  ■  ■  '  ,  ;  t  T  « .  -  '  -  ■  "  /'  1 

r  *  I  _  r  r  y  a  ,  * ■  _  (  n  r  7  a  _  1  -  r  )  -f  S  (  I  9  J  ) 

IF  (  T  •  T  •  ’’I  .OP.  I  •  ^  T  . ' 1 )  CO  TO  14 
I F ( I . FO*  ' I )  GO  TO  2 5 
c ( "?  ,  J ) =-SN 
00  TO  1  4 
?r  S(’"!>J5=-SN 

1  4  50  =  D>rjS  (  GN-S  (  I  9  J  )  ) 

0  T  =  P  7  a  x  1  (  01  ,  SD  ) 

h ' A  —  n p. n S (ON)  v 

RFF  =  DK  4X1  ( REF » HA  ) 

1?  ■ 

c;Frr>  =  ^  j  /ppp 

IF  (  FEP.LE  .ERR.'OR.SK.GP.TIMI  >  GO  TO  IF 
S  K  =  S  K  +  1 

n n  77  |  =  7  , \ ■ 

5  (  j  ,  7 ) =7Pn-2*S (  I  9  p  )  -  F  (  T  9 4  )  /  0  DO  +  S  (  I  *  0  )  /  1  A  00 

2  7  r  (  j  ,  |\  )  =  7  c,  -  ?  0  (  I  » N  I  )  “  S  (  I  9  N  J  )  /  D  0  +  c  (  I  9  N'  <  )  /•!  A  00 

r^0  °P  J=*,NT 

79  r  (  7  ,  J  )  =  "7  P  D  -  7  #  S.  (  G  9 J  )  -  F  (  4  » vT )  /  ?  n  0 + S  (  r  *  J  )  /  1  6  n  0 

7)  =  (_’^r''  -  7^4. 5  (  7  9  '>  )  +S  (  ^  »  7  )  )  -  (  S  (  7  9  4  )  -HF  (  4  9  7  )  1  /  *  o  0  + 

1 ( 5 ( ?  » F ) +0 ( F  ,  7 )  > / 16D0 > /7D0 

f,  (  7  ,  * : )  =  (  7  n  -  ?  *  {  S  (  0  9  M  )  +  S  (  ?  »  M  [  )  )  -  (  S  (  4  9  N  )  +  5  (  7  » N  J  ) ' )  /  ~ 1  0  4 
1  (  S  (  5  >N  )  +  S  (  2  »NK  )) /16D0  )  /2D-0 

GO  TO  1 A 

if  font i nuc 

-J9  IMG  GTREAMRINCTTON  T 0  01  MO  1 1  \/ 

■7  y  =  °po 

r>  f-  p  _  7  r\  7 

17  1*2  *0,1 
II  =  t  +1 
12=1+2 
11=1-1 
T  J  =  I  -  7 
^0  17  J=?»N 
J1=J+1 
J^=J+7 
JI =J-1 
J J= J-2 

TF (  I  .FQ.7 )  GO  TO  18 

VM=  (  5,  (  I  2  9  J  )  -oryo^S  (  II  9  J)  +FOO*S  (  1 1  9  J  )  -S  (  I  J  >  J  )  )  /  (  1  2 ,r'  0  *  H  ' 

GO  JO  10 

;  •  r  (  ^  9  J  )  - 1  9  0 0 *  0  (  -  9  J  )  +  6 HO  S  (  4  *  JO  -S  (  F  9  J  )  )  /  (  1  2  DO#H  ) 
lo  IF ( J  •  FQ  •  2 )  GO  TO  2  0 

Ic (J.FO.R  )  GO  TO  21  . 

l!M=  (  c  (  j  ,  jj  )  (  I  9  J  I  )  +RDO*S  (  I  9  Jl  ) -G  (.1  9  J2  )  )  /  (  1  2f  0*L  > 

GO  TO  22 

j  ;m  _ ( 5 ( t  ,5 ) -fiDO*S ( T  94 ) +18D0*S ( I >3 ) ~1 0D0*$( I  9  7  >  ) / (  1 7  0  0*L  ) 

GC  TO  22  ' 

2  1  ’  !•' 1  ~  (  1  0  0  0  r  (  I  9  r-'  )  ~  1  P  0  0*  S  (  I  9  r '  I  )  +  AO  0 1!-  S  (  I  9  M  J  )  -  0  (  I  9  i*  ))/("'’  ' 1  7C  L  > 
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??  UD=0 An S  (  ■  IM-;  I  (  J  » J  )  ) 

Vn  =  nAR.S  (  V M — V  (  I  ,  J  )  ) 

I =DMAX 1  (  Q I 4 1 JD  4  VO  ) 

HZ'  =n  a  <">s  (  iji^j  ) 

HI'  =  DADS  (  VN  ) 

RFF=DMAX1 (RFF»HA*HR ) 
t:  (  I  *  J)  *1  IN 
17  V (  I  ,  J ) =VN 
VFP=DI/PrF 
OP  TTP  (  A  ,  Ai  )  <  ,  VFP 

I F ( VFR.GF .VTFST . AND • V^R . GE . 5 F-? . ANT . K . GF . 1 n )  CALL  rX I T 

VT  E  c  T  =  VF  7 
FP 1 =rPR*?DO 

IF ( VER.LP.FR 1.0R.K.GF.TI M3 )  GO  TO  2? 

F=K4-1 
GO  TO  ^ 4 

Tyr  PN!D  OF  T  T  rR  A  T  I  ON  S 
7°  -  ’P  T  T F  (  6  ♦  50  ) 


U! R  I  TF  (  A 

*51) 

W  K 

,'4'PR 

'•'R  T  T  r  (  A 

4  A  2  ) 

(  (  W  (  1  , 

J) 

,  T 

—  P 

9  y  i  : 

>  4  J  — 

7 

4  N  ) 

'*'R  T  T  c  (  A 

.  A  7  ) 

WR  ITF (6 

4  51  ) 

r< 

,  yrp 

WRITEi 6 

4  A?  ) 

(  (  T  (  I  * 

J) 

4  I 

=  ? 

, mi  : 

)  4  J  “ 

,!n  ) 

WR T  TF ( A 

,55) 

T  Tr  (  A 

»  M  ) 

SK 

,  S  F  p 

R  I TF ( 6 

4  A?  ) 

(  ( 

S  (  T  , 

J) 

4  I 

=  P. 

, mi  : 

>  ,J  = 

7 

4  N  ) 

K'R  I  IF  (  A 

,  A  4  ) 

'■'?  I  T  c  (  A 

4  A  2  ) 

(  ( 

0(1, 

J) 

4  I 

-P 

, mi  : 

)  »J  = 

7 

.  N  ) 

WR  I  TP ( A 

,  A  7  ) 

(  (  V  (  I  , 

J  ) 

4  I 

_  p 

,  V  1 

>  4  J  “ 

7 

4  M  ) 

1 7  j  Hc  '-'FAN  VELOCITY 
T y  is  THE  mean  T FM0 • 
n T  IS  T Hp  "  T  X F n  vcAN  J F ' • 

A T Ft.  ME  AM,  .  TEMP, MIX  TEMP  4 

STRESS  A  ‘  -  TEMP  GRADIENT  RY  APPLYING  .SIMPSON  RULE 

u 7  •  •  =  •  QC 

T  M  —  r\  Pi  p 

Tw=^no 

r>0  Q  4  J  =  2»N>? 

"  ’  =  no*W  (  y  i  ,  J  ) 

jm-jm  +  t  r>0#T ( 01 , J ) 

? 4  T  '!  =  T  '+?  DO*'--/  ( ■  '  1  *  J  5  *  T  (  M  1  ,  J  ) 

00  7  7  J  =  °  « M  T 9? 

"/  =  K'M+y/  (  y  i  ,  j  ) 

Tm  =  TM  +  T ( '  •  1  . J ) 

7  7  J\  sT’v  +  V  (  0  ,  J)  7-  T  (Ml  ,  J  ) 

OS  7A  I  =  2  4  ' '  4  2 

DO  7  7  J  =  2  » N  » 2 

•••■•■=■■  M+  >■  (  i « j )  *  °  n  o 
t*'  =  tm+t  ( i » j )*?no 
° 7  TV  =  TW+ '  (  T  f  J ) *T ( I  4 J ) *PDP 


« 

• 

4  3 


'  o  6  J-  ■<  ,  M  f  ,  ? 

'*'M  (  I  .  j  )  *  4  r>  o 

Tm  =  Tm  +  t  (  j  .  j  )*4f)n 
^6  TW=TW+H (  T  o J  )#T (  I  * J ) *400 

">0  3  R  T  =3  ,M  J  ,  7 

no  -3  9  J=7,N,3 
|'^=WW+W  (  J  ,  J  )  *4P0 
TM'aTW  +  T  (  I  ,  J  )  *400 
30  TW=TW+W (I»J>*T(I»J) *400 
00  3  8  J  =  ?  T  9 ? 

WM  =  WM+W (  I  9 J ) *?D0 
Tm  =  t'm+t  (  i ,  j  )  *300 
3  8  TW  =  TW  +  W (  I  9 J ) *T (  I  . J ) *3  00 
um_\,/M*?00/  (  ono*F"  1*FM  ) 

TM=TM*?DO/ ( 9D0*FM*FM ) 

TH  =  tv* 3  On/  (  oon*FM*FM  ) 

HR  T  TP ( 6 . 36 ) 

DO  3  0  0  1=1  >  "m 

oh  (  T  ,  1  )  =  (  4  00*  W  (  I  9  3  )  -3  00*0  (  I  )  +4D0*V  (  T. »  4  )  /  3  0  0  ■  (  I  »  8  )  /  4  ; .  o  )  /  L 

TG (  T  9 1  )  =  (  400*  T  (  T  .3  )  —  ^  O  0  *  T  (  T  9  3  )  +40  3*  T  (  I  ,  4  )  /3IOO  -  f  (  T  /41 : 0  )  /L 

OH  (  I  ,  m  1  )  =  (  4-Q0*W  (  T  9 N  )  - 3 i)0* w  (  I  *  N  T  )  +4 00*1'/  (  I  9  N  J  )  /  3 0 0 - 

1  ’■  (  I  ,NK  >  / 40  0  )  /l 

TO  (  T  ,N1  )  =  (400*T  (  I  *N  )  -?DO*T  (hNI)  +40  0*  T  {’i  ,NJ  )  /3D0- 
1 T (  T 9  NK ) / 400 ) /L 
800  roMj  I  m:  if 

00  5  01  J  =  ]  ,  N 1 

•0H(  1  9  J  ).=  (  400* W  (  ?  9  J  )  -3  00*W  {  3  9  J  )  +40 0*0  (  4,  ,  J  )  /  3  0 O  —  (  3  ,  J  )  /4D0  '-/hi 

TG { 1  9  J )  =  ( 400* T ( ? • J ) -3  00*T ( 3 , J ) +400*T ( 4  9 J ) /oOO-T ( 3 , J ) /4O0 >  /h 
3oi  C  O  M  T I M 1 1 F 

103  00  1^4  1=1  9  Ml 

T  TF  (  6 , 3  7  )  T  9  OH  (  T  9  i)  9.3H  (  I  9 N 1  )  9  TG  f  T  9 1  )  9  TG  (  T  ,N1  ) 

104  CONTINUE 

'*'9  T  TP  (  6 , 38  ) 

00  103  J  =  1  ,N1 

WP  I  TF  (f  »  39  )  J  9  .OH  (  1  ,  J  )  ,  TG  (  1  9  J  ) 
loi;  CONTINUE 

OHM  =  SH ( Ml  9 1  ) +SH (Ml  1  ) 

TGM  =  TG ( Ml  9 1 ) +TG (Ml  ,N1  ) 

c,  I J  M  =  0  0  0 

Cl  If  =  Apo 

00  O  o  f  -  3  ,  |V1  ,  p 

SUM  =  SUM+sSH  (  I  9 1  )  +SH  (  I  9  N  1  ) 

0l)T  =  SUT  +  TG(  T  9  1  )  +TG  (  I  *  M  3  ) 
o  o  G ON T  I  N I. ! F 

0  u  h  =  5  H  M  +  4  0  0  *  5  U  M 
T  G  M  =  T  G  M + 4  DO  *  .8 1  )T 

0 1  1  f\/  =000 

OUT  =  00 

00  U.  J  =  3  9M  I  ,  7 

.SIJM  =  SUM  +  SH  (  I  9  1  )  +  S H  (  I  9  N  ]  ) 


c-  U  7  =  C  1  1  '  4-  '  G  (  T  .1  )  +  [  O  (  T  ,  |\|  ]  J 

Q  1.  fON'T  T  f\|  I  IF 

SHMt ( SH^+?DO*SUM ) *H 
TGM=  (  TGm+?D0*SUT  )  *H 
c  ij  m  =  o  D  0 

C  i  I  t  =  n  n  n 

4 

no  ?  ?  j  =  p « m , ? 

SUM=SUM+SH( ] »j) 

SI I T  =  S U T  +  TG  (  1  ,  J  ) 
p  ?  C  O  M  T  I  M 1 1 F 

S  H  M  =  S  H  M  +  su  t  '-*40  n  *  L 
7  G  M  =  T  Gw  +  S 1 1 T  *4  5  0*  L 
Cl  )M=nf)0 

c  i  i  y  =  n  [)  o 

HO  GG  J  =  G  ,N  J  ,  % 

S U  M  =  S  U  M  +  S  H  (  1  ,  J  ) 
ci  '7  =  SHT  +  TG  (  1  i  J  ) 

GQ  COWTIMIIF 

SN»=  (SHU+?DQ*L*Sl  im  )  *?D0*GA  m.va  /  (  300*  (  G  A  "  v.A  +  i  no  )  *'■*?  ) 

TGV=  (  TGM+?r>0*L*SUT  )*?D0*GA  '  -'A/  (  sno*  (G  vima+1  00  )  **?  ) 

FRF  IS  THE  PRODUCT  OF  FRICTION1  FACTOR  AND  REYNOLDS  NO. 
MU  TG  THF  N  ’  l  s'4f  L  T  MO. 

FRF=?QO*SHM/WM 
jn  =  PR*T 

Ml  !  =  V! M *  J  G M  /  T  ’  ■' 

'•’P  T  Tc  (  A  ,  6^  )  ShfM  ,  7GV|  ,  UM  ,  TM  ,  7W 
WPTTf(6,61  )  PRF.TQ.Mi' 

GO  TO  PPG 

4P  FORMAT (  GX  ,  I  ?  »?X  ,  I  G  ,2131  0.  G  ,G  I  p  ) 

46  FORM AT (  PX  ,  D?P . 1 6 , RO] n . Q  ) 

40  FORM A 7 (  1  MJ  ,  GX , ?HM= ,  J  R  ,  GX  ,?HM  =  ,  I ? , GX  ,  PHRFTA  =  .0] a. p , , 

1 QHFRR  ,  n 1  o,3 ,PX , 

1 6 H G A v m A  =  ,D1^.4  ,  P,X  ,3HPR  =  ,  DIO  *  G  ,PX  ,4hRa,C  =  ,D1  0 .3  ) 

^  FORMAT  (  1  H  J  ,  1  0  X  ,  1  7  H  M  0  ^  P  M  T  i  JM  FQlJATKJN) 

p]  FO  R  f1’'  A  T  (  1  OX  ,  6H.SWFFP  =  ,  I  G  0  X  »  4  H  M  A  X  ERR=  ,D13.6  ) 
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Fig.  1  Coordinate  System  and  Finite  Difference  Network  For 

A  Rectangular  Channel 


Fiq .  2  Effect  of  Grid  Size  on  Flow  and  Heat  Transfer  Results  for  Square 

3  4 

Channel  y  =  1 ,  With  Pr  =  0.73  and  ReRa  =  0  and  1.386  x  10 
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Fig.  3  Dimensionless  Axial  Velocity  Distribution  at  X  =  a/2  With  ReRa 
as  a  Parameter  in  a  Square  Channel  y  =  1  and  Pr  =  0.73 
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Fig.  4  Dimensionless  Temperature  Distribution  at  X  =  a/2  With  ReRa 
as  a  Parameter  in  a  Square  Channel  y  =  1  and  Pr  =  0.73 
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Fig.  5  Dimensionless  Axial  Velocity  and  Temperature  Distributions  at  Y  =  b/2  With  ReRa 

as  a  Parameter  in  a  Square  Channel  y  =  1  and  Pr  =  0.73 
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Fig.  6  Streamlines  and  Isothermals  for  a  Square  Channel 
y  =  1  With  Pr  =  0.73  and  ReRa  =  1.030  x  105 
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Fig.  7  Dimensionless  Axial  Velocity  Distribution  at  X  =  a/2  With  ReRa 
as  a  Parameter  in  a  Rectangular  Channel  y  =  2  and  Pr  =  0.73 
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Fig.  8  Dimensionless  Temperature  Distribution  at  X  =  a/2  With  ReRa 
as  a  Parameter  in  a  Rectangular  Channel  y  =  2  and  Pr  =  0.73 
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Fig.  10  Streamlines  and  Isothermals  for  a  Rectangular  Channel  y  =  2  With  Pr  =  0.73  and  ReRa  =  1.408 
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Fig.  11  Dimensionless  Axial  Velocity  Distribution  at  X  =  a/2  With  ReRa 
as  a  Parameter  in  a  Rectangular  Channel  y  =  5  and  Pr  =  0.73 
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Fig.  12  Dimensionless  Temperature  Distribution  at  X  =  a/2  With 
ReRa  as  a  Parameter  in  a  Rectangular  Channel  y  =  5  and  Pr  =  0.73 
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Fig.  13  Dimensionless  Axial  Velocity  and  Temperature  Distributions  at  Y  =  b/2 
With  ReRa  as  a  Parameter  in  a  Rectangular  Channel  y  =  5  and  Pr  =  0.73 
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Streamlines  and  Isothermals  for  a  Rectangular  Channel  y  =  5  With  Pr  =  0.73  and  ReRa  =  2.984 
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Fig.  15  Dimensionless  Axial  Velocity  Distribution  at  X  -  a/2  With  ReRa 
as  a  Parameter  in  a  Rectangular  Channel  y  =  0.5  and  Pr  =  0.73 
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Fig.  16  Dimensionless  Temperature  Distribution  at  X  =  a/2  With  ReRa 
as  a  Parameter  in  a  Rectangular  Channel  y  =  0.5  and  Pr  =  0.73 
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Fig.  17  Dimensionless  Axial  Velocity  and  Temperature  Distributions  at  Y  =  b/2 
With  ReRa  as  a  Parameter  in  a  Rectangular  Channel  y  =  0.5  and  Pr  =  0.73 
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Fig.  18  Streamlines  and  Isothermals  for  a  Rectangular 
Channel  y  =  0.5  With  Pr  =  0.73  and  ReRa  =  1.38  x  10^ 
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Fig.  19  Dimensionless  Axial  Velocity  Distribution  at  X  =  a/2  with 
ReRa  as  a  Parameter  in  a  Rectangular  Channel  y  =  0.2  and  Pr  =  0.73 
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Fig.  20  Dimensionless  Temperature  Distribution  at  X  =  a/2  With  ReRa 
As  a  Parameter  in  a  Rectangular  Channel  y  =  0.2  and  Pr  =  0.73 
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Fig.  21  Dimensionless  Axial  Velocity  and  Temperature  Distributions  at  Y  =  b/2 
With  ReRa  as  a  Parameter  in  a  Rectnagular  Channel  y  =  0.2  and  Pr  =  0.73 
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Fig.  22  Streamlines  and  Isothermals  for  a  Rectangular  Channel 
y  =  0.2  With  Pr  =  0.73  and  ReRa  =  2.535  x  105 
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Fig.  23  w/wn  Versus  ReRa  With  Aspect  Ratio  y  as  a  Parameter 
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Fig.  24  f  Re/(fRe)n  With  Aspect  Ratio  y  as  a  Parameter 
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Table  1  Numerical  Results  for  Square  Channel  y  =  1  and  Pr  =  0.73 
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Table  2  Numerical  Results  for  Square  Channel  y  =  1  and  Pr 
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Table  3  Numerical  Results  for  Aspect  Ratio  y  =  2  and  Pr  =  0.73 
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Table  4  Numerical  Results  for  Aspect  Ratio  y  =  5  and  Pr  =  0.73 
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*■» 

Table  5  Numerical  Results  for  Aspect  Ratio  y  =  0.5  and  Pr  =  0.73 
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Table  6  Numerical  Results  for  Aspect  Ratio  y  =  0.2  and  Pr  =  0.73 
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Table  7  Comparison  of  Nusselt  Number  with  Exact  Value  [20]  for  ReRa 
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Note:  The  above  values  represent  the  greatest  deviation  from  equation  (38). 
For  the  case  y  =  1 »  Pr  =  7.2,  the  differences  are  less  than  the  case 
with  Pr  =  0.73  for  the  same  aspect  ratio. 


